1 Structure

The basic structure of a metric_graph object is defined as

metric_graph <-  R6Class(classname = "metric_graph", 
                         public = list(),
                         private = list(line_to_vertex = function(args),
                                        PtE_to_mesh = function(args),
                                        plot_2d = function(args),
                                        plot_3d = function(args),
                                        coordinates_multiple_snaps = function(args),
                                        merge_close_vertices = function(args),
                                        remove_circles = function(args),
                                        merge.all.deg2 = function(),
                                        remove.first.deg2 = function(args),
                                        compute_lengths = function(args),
                                        A = function(args),
                                        get_edge_weights_internal = function(args),
                                        clear_initial_info = function(),
                                        split_edge = function(args),
                                        compute_laplacian_PtE = function(args),
                                        find_edge_edge_points = function(arga),
                                        add_vertices = function(args),
                                        compute_degrees = function(),
                                        create_update_vertices = function(),
                                        mesh_merge_deg2 = function(),
                                        mesh_merge_outs = function(),
                                        move_V_first = function(),
                                        find_mesh_bc = function(),
                                        set_petrov_matrices = function(),
                                        set_first_weights = function(args)
                                        )
                         )

It only considers classname, public, and private arguments and does not consider any of the other available arguments such as active or inherit. See below.

R6Class(
  classname = NULL,
  public = list(),
  private = NULL,
  active = NULL,
  inherit = NULL,
  lock_objects = TRUE,
  class = TRUE,
  portable = TRUE,
  lock_class = FALSE,
  cloneable = TRUE,
  parent_env = parent.frame(),
  lock
)

2 Some nuances

  • It is not a coincidence that the classname parameter coincides with the variable name of the output of R6Class().

  • The $new() method when we do metric_graph$new() is not defined by the authors but a default one. The $new() method creates the object and calls the initialize() method, if it exists.

3 Unknown (possibly private) functions

  • process_factor_unit()

  • private$set_first_weights()

  • private$line_to_vertex()

  • private$find_edge_edge_points()

4 initialize() inputs

initialize = function(edges = NULL,
                      V = NULL,
                      E = NULL,
                      vertex_unit = NULL,
                      length_unit = vertex_unit,
                      edge_weights = 1,
                      kirchhoff_weights = NULL,
                      longlat = FALSE,
                      crs = NULL,
                      proj4string = NULL,
                      which_longlat = "sp",
                      project = FALSE,
                      project_data = FALSE,
                      which_projection = "Winkel tripel",
                      tolerance = list(vertex_vertex = 1e-3,
                                       vertex_edge = 1e-3,
                                       edge_edge = 0),
                      check_connected = TRUE,
                      remove_deg2 = FALSE,
                      merge_close_vertices = TRUE,
                      factor_merge_close_vertices = 1,
                      remove_circles = TRUE,
                      verbose = 1,
                      lines = deprecated()) {
  # function body
}

5 lines argument


lifecycle::is_present() and lifecycle::deprecate_warn()


if (lifecycle::is_present(lines)) {
         if (is.null(edges)) {
           lifecycle::deprecate_warn("1.2.0", "metric_graph$new(lines)", "metric_graph$new(edges)",
             details = c("`lines` was provided but not `edges`. Setting `edges <- lines`.")
           )
           edges <- lines
         } else {
           lifecycle::deprecate_warn("1.2.0", "metric_graph$new(lines)", "metric_graph$new(edges)",
             details = c("Both `edges` and `lines` were provided. Only `edges` will be considered.")
           )
         }
         lines <- NULL
       }
  • Function lifecycle::is_present() checks if lines was provided by user.

  • Function lifecycle::deprecate_warn() gives a nice formatted warning.

foobar_adder <- function(foo, bar, baz = deprecated()) {
  # Check if user has supplied `baz` instead of `bar`
  if (lifecycle::is_present(baz)) {

    # Signal the deprecation to the user
    deprecate_warn("1.0.0", "foo::bar_adder(baz = )", "foo::bar_adder(bar = )")

    # Deal with the deprecated argument for compatibility
    bar <- baz
  }

  foo + bar
}

foobar_adder(1, 2)
#> [1] 3
foobar_adder(1, baz = 2)
#> Warning: The `baz` argument of `bar_adder()` is deprecated as of foo 1.0.0.
#> ℹ Please use the `bar` argument instead.
#> [1] 3
data_frame <- function(...) {
  lifecycle::deprecate_warn("1.1.0", "data_frame()", "tibble()")
  tibble::tibble(...)
}
df1 <- data_frame(x = 1, y = 2)
## Warning: `data_frame()` was deprecated in <NA> 1.1.0.
## ℹ Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

6 kirchhoff_weights argument


stop()


Here we start to define private fields, which allows to create fields and methods that are only available from within the class, not outside of it. Visit here and here for more information.

if(!is.null(kirchhoff_weights)){
        if(length(kirchhoff_weights)>1){
          warning("Only the first entry of 'kirchhoff_weights' was used.")
          kirchhoff_weights <- kirchhoff_weights[[1]]
        }
        if(!is.numeric(kirchhoff_weights)){
          if(!is.character(kirchhoff_weights)){
            stop("'kirchhoff_weights' must be either a number of a string.")
          }
          if(!(kirchhoff_weights%in%colnames(edge_weights))){
            stop(paste(kirchhoff_weights, "is not a column of 'edge_weights'!"))
          }
        } else{
          if(!is.data.frame(edge_weights)){
            if(kirchhoff_weights != 1){
              stop("Since 'edge_weights' is not a data.frame, 'kirchhoff_weights' must be either NULL or 1.")
            }
          } else{
            if(kirchhoff_weights %%1 != 0){
              stop("'kirchhoff_weights' must be an integer.")
            }
            if((kirchhoff_weights < 1) || (kirchhoff_weights > ncol(edge_weights))){
              stop("'kirchhoff_weights' must be a positive integer number smaller or equal to the number of columns of 'edge_weights'.")
            }
          }
        }
        private$kirchhoff_weights <- kirchhoff_weights
       }
  • Function stop() stops the execution of the entire function at the point where it is located.
sqrt_non_negative <- function(x) {
  print(10)
  if (x < 0) {
    stop("Input must be non-negative.")
  }
  print(25)
  return(sqrt(x))
}

sqrt_non_negative(4)   
## [1] 10
## [1] 25
## [1] 2
sqrt_non_negative(-1)  
## [1] 10
## Error in sqrt_non_negative(-1): Input must be non-negative.

7 tolerance argument

if (is.null(tolerance$vertex_edge) && !is.null(tolerance$vertex_line)) {
           lifecycle::deprecate_warn("1.2.0", "metric_graph$new(tolerance = 'must contain either vertex_vertex, vertex_edge or edge_edge')",
             details = c("`tolerance$vertex_line` was provided but not `tolerance$vertex_edge`. Setting `tolerance$vertex_edge <- tolerance$vertex_line`.")
           )
           tolerance$vertex_edge <- tolerance$vertex_line
           tolerance$vertex_line <- NULL
         } else if(!is.null(tolerance$vertex_edge) && !is.null(tolerance$vertex_line)) {
           lifecycle::deprecate_warn("1.2.0","metric_graph$new(tolerance = 'must contain either vertex_vertex, vertex_edge or edge_edge')",
             details = c("Both `tolerance$vertex_edge` and `tolerance$vertex_line` were provided. Only `tolerance$vertex_edge` will be considered.")
           )
            tolerance$vertex_line <- NULL
         }

        if (is.null(tolerance$edge_edge) && !is.null(tolerance$line_line)) {
           lifecycle::deprecate_warn("1.2.0", "metric_graph$new(tolerance = 'must contain either vertex_vertex, vertex_edge or edge_edge')",
             details = c("`tolerance$line_line` was provided but not `tolerance$edge_edge`. Setting `tolerance$edge_edge <- tolerance$line_line`.")
           )
           tolerance$edge_edge <- tolerance$line_line
           tolerance$line_line <- NULL
         } else if(!is.null(tolerance$edge_edge) && !is.null(tolerance$line_line)) {
           lifecycle::deprecate_warn("1.2.0","metric_graph$new(tolerance = 'must contain either vertex_vertex, vertex_edge or edge_edge')",
             details = c("Both `tolerance$edge_edge` and `tolerance$line_line` were provided. Only `tolerance$edge_edge` will be considered.")
           )
          tolerance$line_line <- NULL
         }

8 which_longlat argument

if(!(which_longlat %in% c("sp", "sf"))){
        stop("The options for 'which_longlat' are 'sp' and 'sf'!")
      }

9 longlat argument

if(longlat){
        private$longlat <- TRUE
        private$which_longlat <- which_longlat
      }

10 proj4string argument

if(!is.null(proj4string)){
        if(!longlat){
          warning("proj4string was passed, so setting longlat to TRUE")
          longlat <- TRUE
          private$longlat <- TRUE
          private$which_longlat <- which_longlat
        }
        private$crs <- sf::st_crs(proj4string)
        private$proj4string <- proj4string
        crs <- private$crs
        private$transform <- !(sf::st_is_longlat(private$crs))
      }

11 crs argument

if(!is.null(crs)){
        if(!longlat){
          warning("crs was passed, so setting longlat to TRUE")
          longlat <- TRUE
          private$longlat <- TRUE
          private$which_longlat <- which_longlat          
        }        
        private$crs <- sf::st_crs(crs)
        private$proj4string <- sp::CRS(crs$input)
        proj4string <- private$proj4string
        private$transform <- !(sf::st_is_longlat(private$crs))        
      }

      if(longlat && (which_longlat == "sp") && is.null(proj4string)){
        proj4string <- sp::CRS("+proj=longlat +datum=WGS84")
        private$crs <- sf::st_crs(proj4string)
        private$proj4string <- proj4string
        private$transform <- !(sf::st_is_longlat(private$crs))        
      }

      if(longlat && (which_longlat == "sf") && is.null(crs)){
        crs <- sf::st_crs(4326)
        private$crs <- crs 
        private$transform <- !(sf::st_is_longlat(private$crs))        
      }

    # private$longlat <- longlat

12 vertex_unit argument

if((is.null(vertex_unit) && !is.null(length_unit)) || (is.null(length_unit) && !is.null(vertex_unit))){
      stop("If one of 'vertex_unit' or 'length_unit' is NULL, then the other must also be NULL.")
    }

    if(!is.null(vertex_unit)){
      vertex_unit <- vertex_unit[[1]]
      if(!is.character(vertex_unit)){
        stop("'vertex_unit' must be a string!")
      }
      if(!(vertex_unit %in% valid_units_vertex)){
        stop(paste("The possible options for 'vertex_unit' are ", toString(valid_units_vertex)))
      }
      private$vertex_unit <- vertex_unit
    }

13 lenght_unit argument

if(!is.null(length_unit)){
      length_unit <- length_unit[[1]]
      if(!is.character(length_unit)){
        stop("'length_unit' must be a string!")
      }
      if(length_unit == "degrees"){
        length_unit <- "km"
      }
      if(!(length_unit %in% valid_units_length)){
        stop(paste("The possible options for 'length_unit' are ", toString(valid_units_length)))
      }
      private$length_unit <- length_unit
    }

if(longlat){
      private$vertex_unit <- "degrees"
      if(!is.null(length_unit)){
        private$length_unit <- length_unit
      } else{
        private$length_unit <- "km"
      }
    } else if(!is.null(vertex_unit)){
        if(private$vertex_unit == "degrees"){
          longlat <- TRUE
          private$longlat <- TRUE
        }
    }

factor_unit <- process_factor_unit(private$vertex_unit, private$length_unit)

14 tolerance argument again

tolerance_default = list(vertex_vertex = 1e-7,
                             vertex_edge = 1e-7,
                              edge_edge = 0)

    for(i in 1:length(tolerance_default)){
      if(!(names(tolerance_default)[i] %in% names(tolerance))){
        tolerance[names(tolerance_default)[i]] <- tolerance_default[i]
      }
    }

    if(verbose > 0){
      message("Starting graph creation...")
      message(paste("LongLat is set to",longlat))
      if(longlat){
        message(paste("The unit for edge lengths is", private$length_unit))
        message(paste0("The current tolerances (in ",private$length_unit,") are:"))
        message(paste("\t Vertex-Vertex", tolerance$vertex_vertex))        
        message(paste("\t Vertex-Edge", tolerance$vertex_edge))           
        message(paste("\t Edge-Edge", tolerance$edge_edge))             
      } else{
        message("The current tolerances are:")
        message(paste("\t Vertex-Vertex", tolerance$vertex_vertex))        
        message(paste("\t Vertex-Edge", tolerance$vertex_edge))           
        message(paste("\t Edge-Edge", tolerance$edge_edge))    
      }
    }    

    if(is.null(tolerance$buffer_edge_edge)){
      tolerance$buffer_edge_edge <- max(tolerance$edge_edge/2 - 1e-10,0)
    }
    max_tol <- max(c(tolerance$vertex_vertex,
                     tolerance$vertex_edge,
                     tolerance$edge_edge))

    private$tolerance <- tolerance

15 V and E arguments


inherits()


Inside methods of the class, self refers to the object. Public members of the object are accessed with self$x, and assignment is done with self$x <- y. Whereas public members are accessed with self, like self$x, private members are accessed with private, like private$x.

PtE_tmp_edge_edge <- NULL
    PtE_tmp_edge_vertex <- NULL

    if(is.null(edges) && is.null(V) && is.null(E)) {
      edges <- logo_lines()
    }
    if(!is.null(edges)){
      if(!is.null(V) || !is.null(E)){
        warning("object initialized from edges, then E and V are ignored")
      }
      if (inherits(edges,"SpatialLinesDataFrame")) {
        tmp_lines = SpatialLines(edges@lines)
        self$edges <- lapply(1:length(tmp_lines), function(i){tmp_lines@lines[[i]]@Lines[[1]]@coords})
      } else if (inherits(edges,"SpatialLines")) {
        self$edges = lapply(1:length(edges), function(i){edges@lines[[i]]@Lines[[1]]@coords})
      } else if(inherits(edges, "MULTILINESTRING")) {
        coords_multilinestring <- sf::st_coordinates(edges)
        lines_ids <- unique(coords_multilinestring[,"L1"])
        self$edges <- lapply(1:length(lines_ids), function(i){coords_multilinestring[coords_multilinestring[,"L1"]==i ,1:2]})
      } else if(is.list(edges)){
        self$edges <- check_lines_input(edges)
      } else {
        stop("edges should either be a list, or of class MULTILINESTRING, SpatialLines or SpatialLinesDataFrame")
      }

    } else {
      if(is.null(V) || is.null(E)){
        stop("You must supply edges or V and E")
      }
      if(ncol(V)!=2 || ncol(E)!=2){
        stop("V and E must have two columns!")
      }
      edges <- list()
      for(i in 1:dim(E)[1]) {
        edges[[i]] <- rbind(V[E[i,1], ], V[E[i,2], ])
      }
      self$edges <- edges
    }
x <- 5
class(x)
## [1] "numeric"
inherits(x, "numeric")  # Returns TRUE, because x is of class "numeric" or inherits from "numeric"
## [1] TRUE
inherits(x, "integer")  # Returns FALSE, because x is not of class "integer"
## [1] FALSE
is(x, "numeric")
## [1] TRUE
is(x, "integer")
## [1] FALSE
is.numeric(x)
## [1] TRUE
is.integer(x)
## [1] FALSE
  • inherits() function is used to check whether an object inherits from a particular class or contains the attributes of a certain class. It returns a logical value (TRUE or FALSE) indicating whether the object inherits from the specified class. It behaves more or less like is() function.

16 set_first_weights() function

self$nE <- length(self$edges)
private$set_first_weights(weights = edge_weights)
  • If everything is ok, set_first_weights() will internally do private$edge_weights <- weights. Press Show below to see the function’s body.
set_first_weights = function(weights = rep(1, self$nE)){
      if(is.null(weights)){
        weights <- rep(1, self$nE)
      }
    if(!is.vector(weights) && !is.data.frame(weights)){
      stop("'weights' must be either a vector or a data.frame!")
    }

    if(is.vector(weights)){
      if ( (length(weights) != 1) && (length(weights) != self$nE)){
        stop(paste0("The length of 'weights' must be either 1 or ", self$nE))
      }
      if(length(weights)==1){
        private$edge_weights <- rep(weights, self$nE)
      } else{
        private$edge_weights <- weights
      }
    } else{
      if(nrow(weights) != self$nE){
        stop("The number of rows of weights must be equal to the number of edges!")
      }
      private$edge_weights <- weights
    }

  }
if(verbose > 0){
      message("Setup edges and merge close vertices")
    }

    t <- system.time(
      private$line_to_vertex(tolerance = tolerance$vertex_vertex,
                           longlat = private$longlat, factor_unit, verbose=verbose,
                           private$crs, private$proj4string, which_longlat, private$length_unit, private$vertex_unit,
                           project, which_projection, project_data)
      )

    if(verbose == 2){
      message(sprintf("time: %.3f s", t[["elapsed"]]))
    }

    if(length(self$edges) > 1){

    if (tolerance$edge_edge > 0) {
    private$addinfo <- TRUE

    if(verbose > 0){
      message("Find edge-edge intersections")
    }

    t <- system.time(
      points_add <- private$find_edge_edge_points(tol = tolerance$edge_edge, verbose=verbose,
      crs=private$crs, proj4string = private$proj4string, longlat=private$longlat, fact = factor_unit, which_longlat = which_longlat)
      )

    if(verbose == 2){
      message(sprintf("time: %.3f s", t[["elapsed"]]))
    }

    PtE <- points_add$PtE

    PtE[,2] <- PtE[,2]/self$edge_lengths[PtE[,1]]

    filter_tol <- ((PtE[,2] > max_tol/self$edge_lengths[PtE[,1]]) &
                     (PtE[,2] < 1- max_tol/self$edge_lengths[PtE[,1]]))

    PtE <- PtE[filter_tol,]

    if(!is.null(PtE)){
      if(nrow(PtE) == 0){
        PtE <- NULL
      }
    }

    if(!is.null(PtE)){
      if(verbose == 2){
        message(sprintf("Add %d new vertices", nrow(PtE)))
      }

      PtE <- na.omit(PtE)

      t <- system.time(
      private$add_vertices(PtE, tolerance = tolerance$edge_edge, verbose = verbose)
      )

      if(verbose == 2){
        message(sprintf("time: %.3f s", t[["elapsed"]]))
      }
    }

    private$clear_initial_info()
    }

    if(tolerance$vertex_edge > 0){
      private$addinfo <- TRUE
      if(verbose > 0){
        message("Snap vertices to close edges")
      }

      t <- system.time(
        PtE_tmp <- private$coordinates_multiple_snaps(XY = self$V,
                                              tolerance = tolerance$vertex_edge, verbose = verbose,
      crs=private$crs, proj4string = private$proj4string, longlat=private$longlat, fact = factor_unit, which_longlat = which_longlat)
        )

      if(verbose == 2){
        message(sprintf("time: %.3f s", t[["elapsed"]]))
      }
      edge_length_filter <- self$edge_lengths[PtE_tmp[,1]]

      filter_tol <- ((PtE_tmp[,2] > max_tol/edge_length_filter) &
                       (PtE_tmp[,2] < 1- max_tol/edge_length_filter))

      PtE_tmp <- PtE_tmp[filter_tol,,drop = FALSE]
      PtE_tmp <- unique(PtE_tmp)
      PtE_tmp <- PtE_tmp[order(PtE_tmp[,1], PtE_tmp[,2]),,drop = FALSE]

      if(!is.null(PtE_tmp)){
        if(nrow(PtE_tmp) == 0){
          PtE_tmp <- NULL
        }
      }

      if(!is.null(PtE_tmp)){
        if(verbose == 2){
          message(sprintf("Add %d new vertices", nrow(PtE_tmp)))
        }

        PtE_tmp <- na.omit(PtE_tmp)

        t <- system.time(
          private$add_vertices(PtE_tmp, tolerance = tolerance$vertex_edge, verbose=verbose)
          )

        if(verbose == 2){
          message(sprintf("time: %.3f s", t[["elapsed"]]))
        }
      }
      private$clear_initial_info()
    }

    if(merge_close_vertices){
      private$merge_close_vertices(factor_merge_close_vertices * tolerance$vertex_vertex, factor_unit)
    }

    if(is.logical(remove_circles)){
      if(remove_circles){
        private$remove_circles(tolerance$vertex_vertex, verbose=verbose,longlat = private$longlat, unit=length_unit, crs=private$crs, proj4string=private$proj4string, which_longlat=which_longlat, vertex_unit=vertex_unit, project_data)
      }
    } else {
        private$remove_circles(remove_circles, verbose=verbose,longlat = private$longlat, unit=length_unit, crs=private$crs, proj4string=private$proj4string, which_longlat=which_longlat, vertex_unit=vertex_unit, project_data)
        remove_circles <- TRUE
    }

    if(merge_close_vertices || remove_circles){
      if(verbose == 2){
        message("Recomputing edge lengths")
      }
      t <- system.time({
        self$edge_lengths <- private$compute_lengths(private$longlat, private$length_unit, private$crs, private$proj4string, private$which_longlat, private$vertex_unit, project_data,private$transform)
      })
       if(verbose == 2){
      message(sprintf("time: %.3f s", t[["elapsed"]]))
       }
    }
    # End of cond of having more than 1 edge
    }

    # Cleaning the edges

    if(verbose == 2){
      message("Post-processing the edges")
    }

    t <- system.time(
          self$edges <- lapply(self$edges, function(edge){
            tmp_edge <- edge[1:(nrow(edge)-1),]
            tmp_edge <- unique(tmp_edge)
            tmp_edge <- rbind(tmp_edge, edge[nrow(edge),,drop=FALSE])
            if(nrow(tmp_edge)>2){
              tmp_edge <- tmp_edge[2:nrow(tmp_edge),]
              tmp_edge <- unique(tmp_edge)
              tmp_edge <- rbind(edge[1,,drop=FALSE], tmp_edge)
            }
            rownames(tmp_edge) <- NULL
            return(tmp_edge)
          }
            )
    )

    if(verbose == 2){
          message(sprintf("time: %.3f s", t[["elapsed"]]))
    }

    # Checking if there is some edge with infinite length
    if(any(!is.finite(self$edge_lengths))){
      warning("There is at least one edge of infinite length. Please, consider redefining the graph.")
    }

    # Checking if there is some edge with zero length
    if(any(self$edge_lengths == 0)){
      warning("There is at least one edge of length zero. Please, consider redefining the graph.")
    }

    end_construction_time <- Sys.time()
    construction_time <- end_construction_time - start_construction_time

    if(verbose > 0){
      message(sprintf('Total construction time: %.2f %s', construction_time, units(construction_time)))

    }

    # Checking if graph is connected
    if (check_connected) {
      g <- graph(edges = c(t(self$E)), directed = FALSE)
      # components <- igraph::clusters(g, mode="weak")
      components <- igraph::components(g, mode="weak")
      nc <- components$no
      if(nc>1){
        message("The graph is disconnected. You can use the function 'graph_components' to obtain the different connected components.")
        private$connected = FALSE
      }
    }
    private$create_update_vertices()

    self$set_edge_weights(weights = private$edge_weights, kirchhoff_weights = private$kirchhoff_weights)

    if (remove_deg2) {
      if (verbose > 0) {
        message("Remove degree 2 vertices")
      }
      t <- system.time(
        self$prune_vertices(verbose = verbose)
      )
      if(verbose == 2){
        message(sprintf("time: %.3f s", t[["elapsed"]]))
      }
    }
    # Adding IDs to edges and setting up their class

    # for(i in 1:length(self$edges)){
    #   attr(self$edges[[i]], "id") <- i
    #   class(self$edges[[i]]) <- "metric_graph_edge"
    # }

    # Cloning the initial graph

    private$initial_graph <- self$clone()

    # Cloning again to add the initial graph to the initial graph
    private$initial_graph <- self$clone()

    # end of function initialize()

17 compute_geodist_PtE() method

compute_geodist_PtE = function(PtE,
                                 normalized = TRUE,
                                 include_vertices = TRUE, verbose = 0){
      if(verbose == 2){
        message("Processing the graph locations...")
      }
      t <- system.time({      
        graph.temp <- self$clone()
        graph.temp$clear_observations()
        df_temp <- data.frame(y = rep(0, dim(PtE)[1]),
                            edge_number = PtE[,1],
                            distance_on_edge = PtE[,2])
        if(sum(duplicated(df_temp))>0){
          warning("Duplicated locations were found when computing geodist. The returned values are given for unique locations.")
          df_temp <- unique(df_temp)
        }

        graph.temp$build_mesh(h = 10000)

        df_temp2 <- data.frame(y = 0, edge_number = graph.temp$mesh$VtE[1:nrow(self$V),1],
                                  distance_on_edge = graph.temp$mesh$VtE[1:nrow(self$V),2])

        df_temp$included <- TRUE
        temp_merge <- merge(df_temp, df_temp2, all = TRUE)

        df_temp$included <- NULL

        df_temp2 <- temp_merge[is.na(temp_merge["included"]),]

        nV_new <- sum(is.na(temp_merge["included"]))

        df_temp2$included <- NULL

        df_temp <- rbind(df_temp2, df_temp)

        df_temp[["__dummy"]] <- 1:nrow(df_temp)

        graph.temp$add_observations(data = df_temp,
                                     normalized = normalized,
                                     verbose=0,
                  suppress_warnings = TRUE)
      })

      if(verbose == 2){
        message(sprintf("time: %.3f s", t[["elapsed"]]))
      }

      if(verbose == 2){
        message("Turning observations of the auxiliary graph to vertices...")
      }      
      t <- system.time(
        graph.temp$observation_to_vertex(mesh_warning = FALSE)
      )
      if(verbose == 2){
        message(sprintf("time: %.3f s", t[["elapsed"]]))
       }
      if(verbose == 2){
        message("Creating auxiliary graph...")
      }
      t <- system.time({
        g <- graph(edges = c(t(graph.temp$E)), directed = FALSE)
        E(g)$weight <- graph.temp$edge_lengths
      })
      if(verbose == 2){
        message(sprintf("time: %.3f s", t[["elapsed"]]))
       }
      if(verbose>0){
        message("Computing geodesic distances...")
      }      
      t <- system.time(
        geodist_temp <- distances(g)
      )
      if(verbose == 2){
        message(sprintf("time: %.3f s", t[["elapsed"]]))
       }      

      if(length(graph.temp$PtV)[1]!=nrow(geodist_temp)){
        un_PtV <- unique(graph.temp$PtV)
        un_coords <- !duplicated(graph.temp$PtV)

        geodist_temp <- geodist_temp[un_PtV, un_PtV]

        tmp_vec <- graph.temp$.__enclos_env__$private$data[["__dummy"]][un_coords]

        un_ord <- order(tmp_vec)

        tmp_vec[un_ord] <- 1:length(tmp_vec)
        #Ordering back in the input order
        geodist_temp[tmp_vec,tmp_vec] <- geodist_temp
      } else{
          geodist_temp <- geodist_temp[graph.temp$PtV, graph.temp$PtV]
          #Ordering back in the input order
          geodist_temp[graph.temp$.__enclos_env__$private$data[["__dummy"]],graph.temp$.__enclos_env__$private$data[["__dummy"]]] <- geodist_temp
      }

      if(!include_vertices){
        geodist_temp <- geodist_temp[(nV_new+1):nrow(geodist_temp), (nV_new+1):nrow(geodist_temp)]
      }

      return(geodist_temp)
  }

18 My change

library(dplyr)
library(here)

# loading the data
load(here("Data_files/data_on_graph_with_covariates_no_consecutive_zeros_partialtomtom.RData"))
load(here("Graph_objects/graph_construction_28_03_2024partialtomtomwhichlonglatsf.RData"))
normalized = TRUE
aux = data_on_graph_with_covariates |>
  rename(distance_on_edge = .distance_on_edge, edge_number = .edge_number) |>
  as.data.frame() |>
  dplyr::select(edge_number, distance_on_edge, .group)

distmatrixlist = list()
PtE = aux %>% filter(.group == as.character(2)) %>% dplyr::select(-.group)


#distmatrixlist[[1]] = sf_graph$compute_geodist_PtE(PtE = PtEE, normalized = TRUE, include_vertices = FALSE)

self = sf_graph
graph.temp <- self$clone()
graph.temp$clear_observations()
df_temp <- data.frame(y = rep(0, dim(PtE)[1]),
                            edge_number = PtE[,1],
                            distance_on_edge = PtE[,2])
if(sum(duplicated(df_temp))>0){
  warning("Duplicated locations were found when computing geodist. The returned values are given for unique locations.")
df_temp <- unique(df_temp)
}
graph.temp$build_mesh(h = 10000)
aa = graph.temp$mesh$VtE
nvertices = sf_graph$nV
#df_temp2 contains all the vertices in graph coordinates
df_temp2 <- data.frame(y = 0, edge_number = graph.temp$mesh$VtE[1:nrow(self$V),1],
                                  distance_on_edge = graph.temp$mesh$VtE[1:nrow(self$V),2])
# creates a new column called included with entries all TRUE
df_temp$included <- TRUE
# just merge. entries in column included from df_temp is TRUE while from df_temp2 is NA
temp_merge <- merge(df_temp, df_temp2, all = TRUE)
# removes column included from df_temp
df_temp$included <- NULL
# this essentially does not change anything
df_temp2 <- temp_merge[is.na(temp_merge["included"]),]
nV_new <- sum(is.na(temp_merge["included"]))
df_temp2$included <- NULL
#df_temp is the same as temp_merge
df_temp <- rbind(df_temp2, df_temp)
df_temp[["__dummy"]] <- 1:nrow(df_temp)
jk = df_temp[,c(2,3)]
nrow(unique(jk))
## [1] 7583

graph.temp$add_observations(data = df_temp,
                                     normalized = normalized,
                                     verbose=0,
                  suppress_warnings = TRUE)
graph.temp$get_data() |> nrow()
## [1] 7583
graph.temp$observation_to_vertex(mesh_warning = FALSE)
graph.temp$nV
## [1] 7581